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SYNOPSIS 

A numerical study of heat transfer from an array of cylinders 
placed in a flow field is presented. Novel techniques such as 
Operator - splitting and Domain Decomposition have been used for 
solving the governing differential equations . Galerkin Finite 
Element Method has been used for accommodating complex geometries 

The main source of error in the solution of heat transfer 
problems occurs in approximating the advective term (U<>VT) at high 
Peclet numbers. In the Operator Splitting algorithm the diffusion 
and advection terms are isolated from the complete equation and 
solved as accurately as possible. Operator Splitting is 
completely free of upwinding and hence free of false diffusion. 
This algorithm is valid for any Peclet number but at high Peclet 
numbers a high order interpolations will give oscillations in the 
temperature field. These have been removed in the present work by 
the Universal Limiters Method. Operator Splitting has a second 
advantage, i.e., it gives reasonable solutions even jon coarse 
grid. This reduces the computational time. 

Practical engineering problems inevitably involve intensive 
computations because of large ^ of physical domains and their 
complex geometries . In many cases computations cannot be 
implemented on existing computers due to limitations of memory. 
The current approach towards handling complex problems is to 
improve the numerical algorithms on one hand and computer 
performance on the other. With the advent of parallel computers 
there exists a great potential in increasing computer performance 
in terms of speed as well as memory. Domain Decomposition 
Techniques are the simplest of parallel algorithms that can 



convert a sequential procedure to one suitable for parallel 
computing. Besides this Domain Decomposition has additional 
advantages which include generation of small matrices and marginal 
reduction in CPU time even on sequential machines. Implicit 
Domain Decomposition and Uzawa's algorithm has been used in the 
present work. 

The model problem taken for analysis is the cooling of an 
array of heated cylinders buried in a porous medium, the earth for 
example. Transient and steady state results have been obtained 
for a Peclet number range of 0.1 to 100. They have been compared 
with analytical solutions whenever possible. 



CHAPTER 1 


INTRODUCTION 

1.1 IMPORTANCE OF PRESENT WORK 

Heat and mass transfer, fluid flow, chemical reactions and 
other related processes occur in engineering equipment, in the 
natural environment, and in living organisms. These processes 
play a vital role in a great variety of practical situations. 
Nearly all methods of power production involve fluid flow and heat 
transfer as essential processes. Major segments of chemical and 
metallurgical industries use components such as furnaces, heat 
exchangers, condensers and reactors, where thermofluid processes 
are at work. Forces that sustain motion of aircraft and rockets 
arise from fluid flow, heat transfer and chemical reaction. In 
the design of electrical machinery and electronic circuits, heat 
transfer is often a limiting factor. The spread of pollution in 
the natural environment is governed by principles of flow, heat 
and mass transfer. 

The prediction of perfomance of a given physical system 
consists of finding the values of the relevant variables such as 
velocity, pressure, temperature and concentration of the chemical 
species in the process of interest at various instants of time. 
The predictions should also describe the effect of changes in 
geometry, flow rates, temperature and fluid properties on the 
performance of the system. Heat transfer and fluid flow processes 
can be studied either from laboratory experiments or by 
theoretical modelling. Experimental studies involve tests 
conducted on a full scale or scale-down model of the physical 
system. However, these tests are in most cases expensive and time 
consuming. Further the small scale models do not always simulate 
all the features of the full scale physical situation. For some 
physical systems, it is not possible to conduct model tests, for 
example, combustion in IC engines. 

A theoretical calculation uses a mathematical model 

consisting of a set of differential equations that describe the 

. . . 

physical system at every point inside the domain and every point 
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in time. Very often, an analytical solution of these equations 
for a practical problem can not be found, thereby limiting its use 
for an engineer. Fortunately, the development of grid based 
numerical methods such as the finite element method and finite 
difference method and the spectacular progress in recent years in 
the speed of digital computers hold promise that the mathematical 
model of almost any practical problem can be worked out. This new 
methodology for attacking the complex problems in fluid mechanics 
and heat transfer has become known as Computational Fluid Dynamics 
and Heat Transfer. The main advantages of such a theoretical 
approach are low cost, remarkable speed with which the numerical 
investigation can be performed, complete information of all 
relevant variables and flexibility to simulate realistic 
conditions of the user's choice. 

Convection-diffusion equations, arising in heat and mass 
transfer is an important class of problems that is being addressed 
in Computational Fluid Dynamics and Heat Transfer (CFDHT) . They 
make appearance in many diverse forms such as Navier- Stokes 
equations with a specified pressure field, vorticity transport 
equation, energy balance eq[uation and chemical transport equation. 
There are many numerical techniques such as finite difference 
method, Petrov-Galerkin finite element method and spectral method 
that are available for solving such equations. 

For more than two-decades the numerical solution of 
convective diffusive transport phenomena has been a very active 
area research for finite difference and finite element modellers. 
The finite element method is very useful for handling complex 
geometries, especially in comparison to the finite difference 
method. The standard Bubnov- Gal erkin finite element approach 
gives rise to spurious oscillations when transport is convection 
dominated. These oscillations result from difficulties associated 
with both the spatial and temporal discretization. The main 
difficulties arising in the numerical solution of the convection 
diffusion equation are due to the non- self adjoint character of 
the associated differential operator. The Galerkin method which is 
well suited to self adjoint problems leads to non-physical spatial 
oscillations when applied to highly convective flows. In such 
problems the convective contribution to the stiffness matrix makes 



it non-syrametric . For getting relatively accurate answer this 
calls for fins' meshes and better (but expensive) matrix inverters. 

For complex geometries finite difference method is difficult 
to use, though the development in finite volume methods 
circumvents the problem to some extent. The finite element method 
is a better choice in this regard. In the finite element context 
most of the efforts to solve convection-diffusion problems have 
used upwind formulations introduced via Petro-Galerkin methods. 
It has now been realized that like their finite difference 
counterparts finite element upwind procedures tend to produce over 
diffusive solutions in multidimensional and transient situations. 
One more disadvantage of upwinding in finite element method (FEM) 
is the problem of selecting proper trial and test functions. 

The above difficulties can be overcome at least partially by 
the operator splitting (OS) method. In the OS formulation, we 
eliminate the convective contribution in Bubnov- Gal erkin finite 
element formulation. This ensures that inversion of highly 
unsymmetric matrix is avoided. The convective term, source and 
decay terms if they are present are solved separately. The method 
of characteristics is used to handle the effect of convection 
alone. Even in the OS method, one gets some undershoots or 
overshoots in the transported variable in highly convective flows 
due to interpolation errors . In the present work the method of 
universal limiters (B. P. Leonard et al, 1990) is used for removing 
inteirpolation related oscillations . 

The present work is concerned with two problems, namely, one 
dimensional unsteady transport and heat transfer from an array of 
heated cylinders placed in a flow field. At high Peclet numbers a 
refined mesh is required over large domains. This necessarily 
leads to the requirement of a large computer memory. The memory 
problem has been resolved in the present study using a domain 
decomposition technique. 

1.2 INTRODUCTION TO OPERATOR SPLITTING 

The operator splitting algorithm (to be called OS) is a 
special case of splitting methods or fractional step methods which 
reduces the solution of a complicated problem to a successive 
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solution of simpler problems. Various kinds of splittings are 
possible. These are given below : 

(a) Geometrical Splitting: This reduces a multidimensional 

problem to a temporal sequence of problems of smaller 
dimensions and in particular to transient one dimensional 
problems. In the latter case we shall speak of splitting in 
terms of spatial direction. A.D.I. (alternating direction 
implicit) scheme for the solution of the transient two 
dimensional heat equation is an example of this kind. 

(b) Physical Splitting: Here the original physical process is 

represented as a temporal sequence of processes possessing 
simpler physical structures. This kind of splitting can~be 
treated as splitting in terras of the physical processes. 

(c) Analytical Splitting: This allows various analytical 

problems to be solved in fractional steps . For example^ the 
restoration of divergence in predictor corrector scheme or 
treating the algebraic terms of the original equations as a 
separate operator in curvilinear coordinates . 

The OS algorithm belongs to the splitting of type (c) . This 
algorithm properly identifies the mixed mathematical character of 
the governing differential ec[uation and solves each homogeneous 
component of the equation as accurately as possible. Hence the OS 
algorithm solves the governing equation true to its character. 

A typical convection-diffusion equation arising in heat 
transfer problems is of the form. 

11 + ( u. 7 T ) = a 7^ T (1.1) 

d d 

In cartesian coordinates V = ( ^ » By ^ velocity vector, T 

is temperature and a is thermal diffusivity. The mathematical 
character of the equation is elliptic if t — > a parabolic if u is 
small and hyperbolic if u is large in magnitude. In all other 
cases, the equation is said to be mixed in character. These 
special cases are summarized below: 


Steady state 

(t — ^ 

> co) (u.VT) = a V^T 

(elliptic) 

(1.2) 

Conduction 

limit 

(lul ■ 

— ^ 0) = a 

(parabolic) 

(1.3) 

Convection 

limit 

{\n\ ■ 

— > co) T^ + U.VT = 0 

(hyperbolic) 

(1.4) 


The OS algorithm treats this problem as follows. At any time 
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level t, one solves Equation (1.4) over a time step At. This is 
considered as a predictor step. The corrector step (Equation 1.3) 
is solved over the same time interval At with the solution 
objtained in the first step as the initial condition. At high 
values of velocity, the correction required is quite small and at 
low velocities, the predictor step acting alone is inadequate. 
However, the two steps confined will furnish the complete solution 
of Equation (1.1) irrespective of the magnitude of u. 

The above approach is completely free of upwinding and hence 
avoids false diffusion errors that usually occur at high Peclet 
numbers. The only error arising in the O.S. algorithm is due to 
time discretization and it is of order At. If either the 
predictor step or the corrector step is solved numerically on a 
grid size of At, then additional truncation errors (~ Ax^, n 2 : 1) 
are possible. Hence the overall error is 

e ~ 0 (Ax^, At) , n ^ 1 

where n depends on the scheme used to solve the diffusive 
corrector step. Numerical techniques are well established for 
solving the diffusion equation. This equation does not have an 
analytical solution when solved in a complex geometry, the finite 
element method is most appropriate for this purpose. In finite 
element method the discretization error depends on interpolation 
function used in the problem. 

The main source of error in solving heat transfer problems 
occurs while approximating the advection terms (u.VT) . In the OS 
algorithm, these terms are isolated from the complete equation and 
solved as accurately as possible. 

On non-dimensionalization, the heat transfer equation becomes 

T^. + (u.VT) = y^T (1.5) 

where Pe = UL/a, is called as Peclet number. It is a measure of 
the relative importance of convective transport with respect to 
diffusive transport. Most of the existing algorithms fail or 
become inaccurate when Pe is raised. The OS algorithm proposed 
here is uniformly valid over any range of Pe. As the Peclet 
Number is raised, a discontinuity in the initial conditions will 
propagate through the flow domain unchanged. This is called a 
front . Predicting the occurrence of a front and its location is a 
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challenging numerical problem. Fronts can also form in non-linear 
problems at low or moderate Peclet numbers. 

The OS algorithm is only conditionally stable since it 
involves explicit time marching and is limited by the Courant 
stability principle. The stability criterion is stated in words 
as, "the numerical movement of the front should not exceed one 
grid spacing in one time increment since the maximum front speed 
is the fluid velocity. Hence, we require ^ ^ - 1/' stability. 

1.3 INTRODUCTION TO DOMAIN DECOMPOSITION 

Practical engineering problems inevitably involve intensive 
computations because of complex geometries and large domains. The 
numerical solution using finite element approximations of 
complicated two and three dimensional nonlinear problems can be a 
formidable task. The difficulties encountered h'kfe are computing 
time and memory space. In many cases the computations cannot be 
implemented on existing computers at all . The current approach 
towards handling complex problems is to improve the numerical 
algorithm on one hand and computer performance on the other hand. 
With the advent of parallel computers there exists a great 
potential in increasing computer speeds as well as memory. 
Parallel computers require parallel algorithm to exploit the 
distributed nature of the processor architecture. Domain 
decomposition techniques are attractive because they can mix 
direct solvers used at the subdomain level and iterative solvers 
used at interface level . Moreover these techniques are ideally 
suited to modern parallel computers, because of the built-in 
parallelism of the algorithm and good localization of the 
associated data. They are fully equivalent to the matrix 
partitioning approach, but have the advantage that subdomains can 
be conveniently chosen based on the problem geometry. Besides the 
ability to parallelize, other specific advantages of domain 
decomposition are given below: 

(1) It can generate smaller matrices that require less memory and 

are rapidly invertible when iterative methods are used. 

(2) A complex assemblage of components of system can be 

systematically analyzed. 



(3) Complex phenomena can be localized by assigning separate 
subdomains to them. 

The basic theme of domain decomposition (DD) is that a large 
domain is divided into many subdomains each of which is easily 
described. The subdomains are linked at the interfaces. The 
governing equations are then solved independently in each 
subdomain. The individual solutions are assembled carefully to 
get a converged solution of the original problem on the complex 
physical domain. The efficiency and accuracy of domain 

decomposition lies in its interface treatment. The subdomains in 
practice may be overlapping or non- overlapping. 

The domain splitting methods originate from well known 
Schwarz method (Le Tallec,1994) for solving elliptic problems. 
First we consider non- overlapping technique with model problem. 

Consider a linear governing differential equation of the type 

Au = f in n (1.6) 


and u = g on 5Q 

where f and g are known functions. Here Q is the physical region 

and an is its boundary as shown in Figure 1.1. A is the Laplacian 

operator. Figure 1.2 shows division of the full domain into two 

subdomains 1 and 2 . The interface between the two subdomains is 

denoted as ' C' . We seek a solution of Equation (1.6) in the 

original large domain. The problem in region 1 is solved 

numerically with an arbitrary Dirichlet boundary condition on C. 

Let this solution be u^, u being the original dependent variable. 

This is followed by numerical solution of u in region 2 with a 

boundary condition on C which is derived from u^ . The global 

solution u and its derivative should be continuous at the 

interface. Hence the most obvious choice of boundary condition on 

1 


C for region 2 would be normal derivative 
the normal at C exterior to region 1. 

solved in region 2 . 

In mathematical form, we have in region 1 


5u^ 

an 


Here n denotes 


With this boundary 
condition u can be solved in region 2. Let this solution be u^. 


Au^ = f in (1.7) 

u., = g on 3Q, and u., = 

1 1 1 


X on C 
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In region 2 



= f in ^2 

g on 8Q^ and 




at C 


Here SQ = an^ U an^ and Q ~ OJ 

Parameter X is the Dirichlet boundary condition for 
described earlier and n^ and n 2 are the normals at C exterior to 
regions 1 and 2 respectively. 

The problem in region 1 is solved again with a new boundary 

condition on C. This may be the value of u^ on C or some 

11 . 2 
combination of U 2 and u^ on C. Let the solution be denoted u^. 

The procedure continues in this fashion where the problem in 

regions 1 and 2 are solved alternately until convergence in the 

global solution is achieved. This convergence criterion will 

satisfy the conditions of continuity of u and its derivative at 

the interface. The method described above is known as the 

alternating Dirichlet -Neumann (D-N) Method because Dirichlet and 

Neumann boundary conditions are applied alternatively on the 

interface . 

The domain decomposition method described above is sequential 

in nature. Its efficiency depends on how the iterations are 

handled at the subdomain interfaces . Consider the alternating D-N 

method as an example. After one complete sweep of the domain, the 

use of arithmetic means of u- ^ and u„^ as the boundary condition 
I+l ^ 

on C for u^ can lead to a considerable improvement in overall 
convergence rate. Here notation I represents iteration number. 
Thus the interfaces are the most crucial zones and they should be 
handled with care and intelligence. In the present work, Uzawa's 
algorithm has been used for interface iterations . ^ 


1.4 LITERATURE SURVEY 

The operator splitting technique, which is being presented 
here as an alternative to the upwind method was first proposed by 
Yanenko et al. (1984) . Various possible types of operator 
splittings feasible in gas dynamics equations have been described 
in their work. Issa (1985) describes a non iterative method for 
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handling the coupling of the implicitly discretized time dependent 
fluid flow equations. The method is based on the use of pressure 
and velocity as dependent variables and hence applicable to both 
compressible and incompressible versions of the momentum 
equations . The main feature of the technique is the splitting of 
the solution process into a series of steps whereby operations on 
pressures are decoupled from those on velocity at each step, with 
split sets of equations being amenable to solution by standard 
techniques. Cooke (1986) has mathematically established that a 
certain operator splitting of the two dimensional Navier-Stokes 
equations in conservative form is second order accurate not only 
for linear systems, but also is true in the presence of 
non-linearity. Glowinksi et al. (1992) describe various operator 
splitting methods for advection - diffusion equations and their 
drawbacks . 

Ding and Liu (1989) present a new operator splitting 
algorithm for the two dimensional convection-dispersion-reaction 
equation. The governing equation is split into three successive 
.initial value problems : a pure convection problem, a pure 
dispersion problem and a pure reaction problem. For the pure 
convection problem solutions are found by a high order upwind 
procedure. For the pure dispersion problem, time-explicit finite 
element algorithm is employed. Analytical solutions are obtained 
for the pure reaction problem. Results presented are mainly for 
one dimensional problems . 

Muralidhar et al . (1993) have applied the operator splitting 
algorithm to two dimensional problems in simple geometries and 
describe the advantages of operator splitting over upwinding. 
They have subsequently applied it to enhanced oil recovery using 
hot water injection (Muralidhar and Pillai, 1993) . 

Leonard and Mokhtari (1990) present a method for removing 
oscillations encountered in high order upwinding schemes. Results 
presented are mainly for one -dimensional problems. Leonard and 
Nikhafs ri991) discuss a strategy for accurately simulating highly 
convective flows containing discontinuities such as density fronts 
or shock waves^ Without distorting the smooth profiles or 
clipping narrows local extreme and disadvantages of other shock 
capturing methods. This method has also been applied only to 
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one -dimensional problems. Westerink and Shea (1989) explain how 
to handle oscillations arising due to strong convection in finite 
element method through the use of a Petrov-Galerkin method. 

Glowinksi et al (1983) discuss domain decomposition methods 
for nonlinear problems in fluid dynamics in complex geometries 
using finite element method. They have given a summary of domain 
decomposition methods suitable for parallelization using Schwarz's 
algorithm. 

Yagawa et al (1991) have developed a parallel finite element 
method based on domain decomposition for stress analysis. They 
use Uzawa's algorithm and the conjugate gradient method for 
interface treatment. Chatter jee (1993) has applied domain 
decomposition technique to oil recovery and transient non-linear 
heat conduction. He uses Uzawa's algorithm for interface 
treatment . 

Laura and Habashi (1994) present the parallelization aspects 
of a solver for fully coupled 3-D compressible Navier — Stokes 
equations. Finite element method is used for solving the 
differential equations and conjugate gradient method is used for 
interface treatment. They explain the difficulties encountered in 
domain decomposition methods owing to book keeping of geometry and 
convergence criterion at interfaces for non-linear problems. 

1.5 SCOPE OF THE PRESENT WORK 

The present work is in two parts. In the first part the 
performance of operator splitting for advection-dif fusion problem 
is tested for simple as well as complex geometries. In the 
present work, the velocity field is calculated using Darcy's law 
applied to flow in a porous medium. With operator splitting, the 
diffusion equation is solved using finite element method and the 
advection step is solved using method of characteristics. The 
results obtained with operator splitting are compared with 
analytical results whereever possible. To remove errors due to 
higher order interpolation, universal limiters have been used. ' 

In the second part of this study, the performance of domain 
decomposition is tested for simple as well as complex geometries. 
Only diffusion step is solved locally in each subdomain, and 
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advection step is solved globally over the whole domain, 
algorithm is used for interface treatment . The effect of 
of outflow plane has been studied using two subdomains, 
present work, the original domain is divided into 2, 
subdomains and their performance has been studied. 


Uzawa' s 
location 
In the 
4 and 9 




CHAPTER 2 


MATHEMATICAL FORMULATION 


2.1 GOVERNING DIFFERENTIAL EQUATIONS, BOUNDARY CONDITIONS AND 
INITIAL CONDITIONS 

Consider the differential equation governing energy 
transport, 

Tj. + u.VT = a V^T (2.1) 

where T is temperature, u is velocity vector and a is thermal 

diffusivity. Non-dimensionalization of Equation (2.1) is carried 

out by using R as the length scale, R /a as time scale, approach 

velocity U as the velocity scale and qR as the temperature scale. 

For flow past a cylinder or an array of cylinders, R is the 

cylinder radius and q is the specified heat fliix on its surface. 

In a problem without cylinders, R is a longitudinal distance and 

the temperature scale is taken as T - T . . 

^ max min 

Equation (2.1) in dimensionless form is, 

T^. + u.VT = V^T (2.2) 

where Pe = UR/a is Peclet number. 

The boundary conditions chosen for the present study are given as 
follows . 

Case X: Movement of a temperature front in a semi infinite medium 

x=0, 0:£y:sL;T = l 

X = L, 0 y L; T = 0 

y = 0 and L, O^xsL, |^=0 

where n is a unit outward normal and L is length scale . 

The initial condition is taken asT = 0att = 0, i.e., the 
medium is initially cold. We study the response of the medium 
when it is subject to hot fluid suddenly at its inflow plane. 

Case II : Flow Past an Array of Heated Cylinders 


X = 0, 
X = L, 


0 ^ Y ^ L; 
0 s y :s L; 
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y = 0 and L, 0 ^ x s L; 55 . 0 

5T 

On the surface of the cylinders, ^ 
where n is a unit outward drawn normal. 

The initial condition is T = 0, at t = 0, as before. The 
velocity field in both problems is steady and prescribed. 


2.2 OPERATOR SPLITTING FORMULATION 

The OS algorithm properly identifies the mixed mathematical 
characters of the governing differential equation and solves each 
equation as accurately as possible. 

Consider the governing differential equation for energy 
transport in 2 -dimensional 


T^ + U.VT = V^T 


(2.3) 


Application of the OS algorithm to equation (2.3) generates the 


following steps . 


1 2 

T = — V T 
t Pe 


(2.4) 


T^ + U.VT = 0 (2.5) 

At any time interval, we solve equation (2.4) over a time 
step At. This is the predictor step. The corrector step is 
subsequently solved over the same time interval At with the 
solution of the predictor step as the initial condition. The two 
steps combined will completely solve the full equation 
irrespective of the magnitude of Peclet number. In the present 
study, the predictor step is solved using a Bubnov-Galerkin finite 
element method. The finite element method formulation is 
described in Appendix A. The corrector step is solved using the 
method of characteristics and is described in Appendix B. 


2.3 UNIVERSAL LIMITERS 

Accurate simulation of highly convective flows continues to 
be one of the most challenging problems in computational 
mechanics. The inability to adequately simulate simple scalar 
profiles even in straight forward case of 1 -dimensional pure 
convection at constant velocity has been referred to as the 
ultimate embarrassment for computational mechanics . All existing 
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methods are well suited for handling smooth solutions . A sudden 
change in the value of dependent variable gives spurious 
unphysical oscillations. 

-Jn the present work no upwinding has been used. However, in 
the operator splitting algorithm a high order interpolation 
procedure is used for calculating the temperature . Interpolation 
methods of order two or higher are themselves a source of 
oscillations. To limit the oscillations, the method of universal 
limiters (Leonard & Niknafs, 1991) has been adopted in the present 
work and is described below. 

Figure (2.1) shows the local behaviour of any convected 
variable (p . Depending on the direction of the local velocity we 
label the three indicated node values (upwind), (downwind) 
and 4>q (centrally located between (p^ and <P-q) • In Figure (2.1(a)) 
the local variation of (p is monotonic. One necessary condition of 
the universal limiter is sketched in this figure and is shown by 
the cross-hatched limits. The convected face value <p^, should lie 
between adjacent node values in locally monotonic regions. 

The treatment of locally non-monotonic behaviour is shown in 
Figure (2.1) (c) and (d) . Handling this case is more difficult 
compared to monotonic behaviour. 

It is convenient to summarize the limiter constraints in 
terms of normalized variables as follows. Let 



Hence = 0, = 1 


Symbolically, the' universal limiter constraints on can be 
written (in monotonic case) as. 


?c - ^f - 1' 0 s i 1 (2.5) 

and. 0 s s 

This above condition is not enough to guarantee computational 
monotonicity. For this consider the explicit update of ? under 
conditions of constant velocity. In terms of normalized 
variables, 

? ^ ^ - C (?^ - ^ ) 

^c ^c ^f ^u 
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where C = — 7 = is the Courant number. 

^ In order to maintain monotonicity (locally) , must satisfy 

^ n +1 ^ ^ n +1 ^ . n+l 
^u ^c 

Taking the conservative worst case estimate implies, 

0 ^ ^ ^ 1 
c 

From Equation (2.5), the right hand side of the inequality is 
assured. The left hand inequality leads to, 

^f ^u ^c' 

Again the worst case estimate gives = 0 snd so, 

^ ( 2 . 6 ) 

Thq inequalities Equation/)(2 . 5) and Equation (2.6) constitute 
the universal limiter constraints on with respect to when 
is within the monotonic range (0 < 1 ). Figure (2.2) shots^a 

clear picture of these constraints . The algorithm for applying 

the method of universal limiters is as given below. 

1 . According to the direction of velocity identify <p , <p and 

\x c 

2. Compute (p^ using some interpolation scheme, say cubic or 
quadratic. 

3. Compute the corresponding normalized variables. 

4. Check the limits on (Equation (2.5), and Equation (2.6)). 

5. If falls within the region (Figure 2.2) skip the limiter's 

algorithm and go to , 8 , if not go to 6 . 

6 . Choose the nearest possible value of ? ^ . 

7 . Denormalize with following expression 

<^f = - <P^) + (P^ 

8 . Repeat this procedure for each node . 

The above approach is given for a one -dimensional problem. 
In the two-dimensional case, the algorithm is implemented along a 
streamline assuming it to be locally straight. 

2.4 FORMULATION FOR DOMAIN DECOMPOSITION 
2.4.1 Uzawa's Algorithm with an Example 

We first outline Uzawa's algorithm for an elliptic problem to 




show its parallel nature and hence provide a clear motivation for 
its use in the present study. Consider a linear second order 
partial differential equation of the type, 

Au = f in n and u=gonSQ (2.7) 

Here f and g are known functions, fi is the physical region and dQ 

is its boundary (Figure 1.1). Figure 1.2 also shows a decomposition 

of Q into two subdomains 1 and 2 with C as their interface. For a 

globally converged solutiqn one requires both u and its normal 

derivatives to be continuous on C. We start with a guess value of 

the derivatives, I and I ( = x^) where n., and n- are the 

3n^ c Sn 2 c 12 

normals at C exterior to regions 1 and 2 respectively. Let this 

guess be . Using X"^ as the boundary condition on C and solves 

the governing differential equation in regions 1 and 2. Let the 

■ 1 1 1 2 
solutions respectively be u^ and u^ . Let y^ and be the values 

11 1 
of u^ and U 2 ^ respectively on C. Then X is updated as, 

^2 .1 ,1 T, - 

X = X + p (y^ - y^) /L 

Here p is an adjustable constant that lies between 0 and 1 and L 

2 

IS a suitable length scale. Using X as the next guess, the 
entire process is repeated. After I iterations, Uzawa's Algorithm 
for updating the interface flux is, 

= X^ + p(y2^ - y^^)/L (2.8) 

Interface convergence is judged by the criterion 

- X^l s 5 

where 5 is the specified convergence griterion. 

The iterations of the algorithm given above converge when the 
computed function u^ (x) matches the steady state solution of 
Equation (2.7). In a transient problem, the same algorithm as 
applied between successive time steps and each subproblem in the 
subdomains are solved by an implicit finite element method, The 
choice of the implicit scheme guarantees the stability of the 
time-marching procedure and convergence is dictated by the 
magnitude of interface parameter, p alone. When the gradient of 
the dependent variable converges, the variable at the interface 
becomes continuous and this is implied by Equation (2.8). Hence 
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the algorithm ensures continuity of u as well as its normal 
derivative at the interface. It is worth noting that iterations 
are completely avoided if the transient problem is solved by an 
explicit time-marching procedure. 

It can be seen that Uzawa's algorithm given above is 
completely parallel in nature. Calculations in regions 1 and 2 
can be carried out independently on two different processors of a 
parallel computer. The processors return the value of y^ and y^ 
after each iteration to b^ie front end processor (FEP) . Tlie FEP 
manipulates the incoming data and performs convergence tests 
before sending the data back to the individual processors for 
subsequent calculations. In general, one can handle as many 
subdomains as the number of processors available. 

Equation (2.8) showB that domain decomposition introduces 
interface iterations even in linear problems . Hence the 
possibility of an increase in CPU time must be considered, 
especially under transient conditions. However results presented 
by Chatter jee (1993) show that CPU time marginally decreases with 
increase in the number of subdomains even on sequential machines. 
The reason for this result is the following. When matrix 
inversion is accomplished using iterative methods or direct 
methods e.g. Gauss-Seidel or Gauss elimination,- the convergence 
rate varies as N , where N is the matrix size and a is some 
positive constant. Hence convergence is significantly delayed for 
large values of N. Application of domain decomposition cuts the 
matrix size by a factor greater than or equal to 2, depending on 
the number of subdomains used. Hence each subdomain calculation 
is accelerated, providing an overall computing time advantage. 

2.4.2 Movement of a Front: 2-subdomains 

1 2 

Governing differential equation: T^ = ^ V T 
0 < X < 18, 0<y<6, t>0 

Initial conditions: T = 0 at t = 0 


Boundary conditions :x=0, 0^y^6;T=l 

X = 18, 0:sy:£6;T = 0 
y = 0 and 6, 0 x ^ 18 ; = 0 


(2.9) 



The region 


As before, n is outward drawn nomnial on a surface, 
size 18 X 6 has been arbitrarily chosen in this work. 

Figure (2.3) shows the decomposition of the original problem 
into two equal subdomains . While the interface can be chosen to be 
an arbitrary surface, it is most convenient to have it overlap 
with one of the coordinate lines. In the present example, it is 
the coordinate line x = 9. The mathematical problem solved in 
each subdomain are given below. 


Subdomain 1 

The governing differential equation: 

0<x<9, 0<y<6, t>0 

Initial condition: T = 0 at t = 0, 



( 2 . 10 ) 


Boundary conditions: x 

X 

y 


0, 0sy:s6;T = 1.0 


e. 0 S y ^ 6; II = - 


aT 


0 and 6, 0^x£9;^=0 

on 


Subdomain 2. 

1 2 

Governing differential equation: T^ = ^ V T, 

6 < X < 12, 0<y<6, t>0 (2.11) 

Initial condition: T = 0 at t = 0, 

3T I 

Boundary conditions: x = 6, 0:sy<6;^=X 

x = 18, 0:£y^6;T = 0 
y . 0 and 6, 9 s x s 18; || = 0 ■ 

Equations (2.10) and (2.11) are solved by an implicit Bubnov 
Galerkin finite element method using six-noded isoparametric 
elements. The stiffness matrix is inverted using a sparse matrix 
solver (MA28) . The matrix solver uses Gaussian elimination 
technique with full pivoting. At the I-th iteration between 
successive time steps, let T^ and T^ be the teirperature values 
computed at the interface from subdomain 1 and 2 respectively. 
Then, the guessed gradient boundary condition X is updated using 
Uzawa's algorithm as follows. 

= a" * P (t/ - T^h/L 


( 2 . 12 ) 
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For a characteristic dimension of L = l, p = 0.2 has been 

found to yield optimum results in the earlier work of Chatterjee 

(1993) . This value has been used again in this study. Figure 

(2.4) shows a flow chart for implementing Uzawa' s algorithm. The 

convergence criteria used is - X^l :£ 5, where 5 is a small 

-4 . 

value, typically 10 in dimensionless calculations. 


2.4.3 Movement of a Front: 4-Subdomains 

. . 1 2 
Governing differential equation: T^ = ^ V T 

0 < X < 18, 0<y<6, t>0 (2.13) 

Initial condition: T = 0 at t = 0 

Boundary conditions: x=0, 0sy<6;T=l 

x=18, 0sy^6;T=0 

y = 0 and 6, 0:sxil8;|^=0 

Figure (2.5) shows the decomposition of the original problem 
into four equal subdomains. The mathematical problems solved in 
each subdomain are given below. 

Subdomain 1. 

. . 1 2 
Governing differential equation: T^ = ^ ^ 

0<x<9, 0<y<3, t>0 (2.14) 

Initial condition: T = 0 at t = 0 



Initial condition: T = 0 at t = 0 

5T I 

Boundary conditions: x = 9, 0^ys3;^=X^ 



Subdomain 3 


x = 18, 0 :sy: 53 ;T = 0 
y = 0, 9 5 X i 18; = 


0 


y = 3, 9 ^ X ^ 18; II = -xl 

1 2 

Governing differential equation: ^ \7 T 

0<x<9, 3<y<6, t>0 

Initial condition: T = 0 at t = 0 

Boundary conditions: x = 0, 3:£y:£6;T=l 

X = 9, 3 ^ y ^ 6 ; ^ = - X 3 


X < 9- ^ 
^ an 


aT 


y = 3, 0 

y = 6 , 0 ^ X ^ 9; 

Subdomain 4 

Governing differential equation: 

9 < X < 18, 3<y<6, t>0 

Initial condition: T = 0 at t = 0 


= 0 


1 2 

T = ^ V^T 
t Pe 


X = 9, 

7 < fi. ST _ 

- y - Q^- 

4 

X = 18, 

3 ^ y ^ 6; T = 

0 

y = 3, 

y - X - la, - 

^ 2 ^ 

y = 6 , 

9 X 18 • — - 

y - X - la, - 

0 


(2.16) 


(2.17) 


At the I-th iteration between successive time steps, let T 


and T 


14' 


■"21 ^ 22 ' 


T34 and T33, 


and T ^3 and T ^2 


11 
be the 


temperature values computed at interfaces from subdomains 1 , 2, 3 
and 4 respectively. Then the guessed gradient boundary condition 
X is updated using Uzawa' s algorithm as follows . 


i+i 

1 


+ 

p 

(T^ 

^21 

- ) 
11^ 

/ 

L 

(2.18) 

i+i 

2 


+ 

p 


- ) 
22^ 

/ 

L 


i+i 

3 

" ^3 

+ 

p 


- ) 
J. 33 ; 

/ 

L 


i+i 

4 

= ^4 

+ 

p 

*'^34 

- ) 
14^ 

/ 

L 



As stated earlier, values of L = 1 and p = 0.2 have been 



used in the present work. The flow chart for implementing Uzawa's 
algorithm is similar to the one given in the context of 2 
subdomains . 


2.4.4 Movement of a Front; 9 - Subdomains 

1 2 

Governing differential equation: ^ V T 

0 < X < 18, 0<y<6, t>0 (2.19) 

Initial condition: T = 0 at t = 0 


Boundary conditions: x=0, 0:£y:s6;T=l 

x = 18, OiS y£6;T=0 

QI-p 

y = 0 and 6, 0:£x:si8;^=0 

Figure (2.6) shows the decomposition of original problem into 
9 equal subdomains . The mathematical problem solved in each 
subdomain are given below. 


Subdomain 1 

Governing differential equation; 


1 2 

T ^ ^ 7 T 
t Pe 



0 < X < 6, 

O 

A 

< 

2 

/ 

t 

> 

0 


Initial 

condition: 

T = 0 

at 

t 

= 

0 




Boundary 

conditions : 

X = 0, 

0 

;< 

y 

< 

2; 

T = 

1 



X = 6, 

0 

:< 

y 

< 

2; 

aT 

an " 

1 

H H 



il 

o 

0 

< 

X 

;< 

6; 

aT 

an " 

0 



y = 2, 

0 

< 

X 

< 

6; 

aT 

an 



Subdomain 2 

Governing differential equation: 

0<x<6, 2<y<4, t>0 

Initial condition: T = 0 at t = 0 


1 2 

T = ^ V T 
t Pe 


il 

O 

to 

:< 

y 

< 

4; 

X = 6, 

2 


y 

< 

4; 

II 

to 

0 


X 

s 

6 ; 

il 

0 


X 

< 

6; 


4; T = 1 

£1 
an " 
aT 


.1 

an “ ^2 


aT 

an 


I 

^4 


( 2 . 20 ) 


( 2 . 21 ) 



Subdomain 3 


1 2 

Governing differential equation: ^ V T 

0<x<6, 4<y<6, t>0 

Initial condition: T = 0 at t = 0 

Boundary conditions: x = 0, 4:sy<6; T=1 

X = 6 4 — y — 6* 
y = 4, 0 = X ^ 6; II » ^4 


0 


1 2 

T = — V T 
t Pe 


y = 6, 0 ^ X ^ 6; II = 

Subdomain 4 

Governing differential equation: 

6<x<12, 4<y<6, t>0 

Initial condition: T = 0 at t = 0 

5T 

Boundary conditions: x = 6, 4sy:£6;^ = A. 


x = 12, 4^y£6;-^=-/\. 


81 

3n 


y = 4, 6 ^ X ^ 12; II = xj 


y = 6, 6 s X ^ 12; 


an 

81 

an 


( 2 . 22 ] 


(2.23) 


Subdomain 5 

. . . 1 2 
Governing differential equation: ^ V T 

12 < X < 18, 4<y<6, t>0 

Initial condition: T = 0 at t = 0 


X = 12, 4 

VI 

VI 

6; 

aT 

an ~ 

x^ 

X = 18, 4 

VI 

VI 

6; 

T = 

0 

y = 4, 12 

VI 

VI 

18 

aT 
' an 

M CO 

II 


y = 6, 12 s X ^ 18; II = 0 


Subdomain 6 


1 2 

Governing differential equation: V T 

12<x<18, 2<y<4, t>0 


(2.24) 


(2.25) 
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Initial condition: 
Boundary conditions : 


T = 

0 

at 

t 

= 

0 




X = 

12 

, 2 


y 

:S 

4; 

ax 

an 

x^ 

X = 

18 

, 2 

< 

y 


4; 

X = 

0 

y = 

2, 

12 

:< 

X 

< 

18 

ax 
' an 

O 
H H 

il 

y = 

4, 

12 

S 

X 

< 

18 

ax 
' an 

H 00 

1 

II 


Subdomain 7 

Governing differential 

12 < X < 18, 

Initial condition: T 

Boundary conditions : x 

X 

y 

y 


. 1 2 
equation: “ Pe ^ 

0<y<2, t>0 

= 0 at t = 0 

= 12 , 0 » y * 2 ; II = 

= 18, 0£y:£2;T = 0 
= 0, 12 X :£ 18; = 0 


2, 12 


18; 


5n 


I 

10 


Subdomain 8 

Governing differential equation: 

6<x<12, 0<y<2, t>0 

Initial condition: T = 0 at t = 0 


1 2 
_ 

Pe 


X = 6, 

Yl 

>1 

VI 

O 

2; 

W X 

X = 12 , 

o 

lA 

^ 2; 

ax 

an “ 

11 

o 

6 ^ X ^ 

12; 

ax 

an “ 

y = 2, 

6 X ^ 

12; 

ax 

an 


I 

^11 


I 

'12 


Subdomain 9 

Governing differential 
6 < X < 12 , 

Initial condition: T 

Boundary conditions : x 


1 2 

equation: “ Pe ^ 

2<y<4, t>0 


= 0 at t = 0 


6, 2 £ y 4; 
12, 2 s Y s 4; 


ar 

an 

ax 

an 


= x: 


(2.26) 


(2.27) 


(2.28) 


X 





Figure 2.5 Division of the Computational Domain into Fi 
Subdomains 
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Figure 2.4 Flow chart for Domain Decomposition with Uzawa' 
Algorithm 
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Figure 2.6 Division of the Computational Domain into Nine 
Subdomains 
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y = 2 , 
y = 4 , 6 ^ X 


5 < X ^ 12 • — = 

b - X - - 


12 ; 


aT 

an 


X 

= - X 


I 

12 

I 

6 



At 

the I- 

-th iteration between 

successive 

time 

steps, let 

T-l 

11 

and 

12 ' 

22 ' 

T ^3 and T 34 and 

35' 

T 45 ' 

^46 “1 ^47' '^57 

and 


tI 

68 

' 69 

■r?io 

■^711' 

rpl 

811' 

812 

and 

npl 

93' 


III 

"^gg' "^99 '^912 temperature values computed at the interfaces 

from subdomains 1, 2, 3, 4, 5, 6 , 7, 8 and 9 respectively. In the 

X represents subdomain number, y represents 


notation 


xy' 


interface number and I represents iteration number. The guessed 


gradient boundary condition X 
as follows 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 

X 


I+l 

T 

= 

^ 1 ^ 
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p 

(T^ 
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/ 

L 
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= 

^2 
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p 
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'22 

) 

12^ 

/ 

L 

I+l 
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^^1 

4 

p 

(T^ 
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(T^ 
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'^812* 
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is updated using Uzawa' s algorithm 


(2.29) 


L 

L 

L 


2.4.5 Heat Transfer from a Single Cylinder Using 2-Subdomains 

. 1 2 

Governing differential equation: “ Pe ^ 

0 < X < 12, 0<y<6, t>0 

Initial condition: T = 0 at t = 0 

Boundary conditions: x=0, 0:£y£6;T=0 

X = 12, 0 ^ y = 6; = 0 


(2.30) 
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dT 


0 


y = 0 and 6, 0 :£ x :s 12; 

On the cylinder surface, we consider a constant heat flux 
boundary condition 
- 1 , 


aT 

an 


where, n is a unit outward drawn normal on the cylinder surface. 

Figure (2.7) shows the decomposition of original problem into 
2 equal subdomains. Each subdomain has a size of 6 x 6 units. The 
mathematical problem solved in each subdomain are given below. 


1 2 

T = — V T 
^t Pe 


Subdomain 1 

Governing differential equation: 

0<x<6, 0<y<6, t>0 

Initial condition: T = 0 at t = 0 

Boundary conditions: x=0, 0£y^6;T=0 

aT 


(2.31) 


X = 6, 0 


6 ; 


- X" 

ar 

an 


- - an 

y = 0 and 6, 0 ^ x ^ 6; 
On the cylinder surface, - ^ = 1. 

Subdomain 2 

1 2 

Governing differential equation: T^ = ^ V T 
6<x<12, 0<y<6, t>0 

Initial condition: T = 0 at t = 0 

5T I 

Boundary conditions: x=6, 0£ys6;T::r = X 


= 0 


(2.32) 


an 

12 0 ^ y ^ 6 • — = 

1 ^, U y an 


aT 


0 


y = 0 and 6 , 6 ^ x ^ 12 ; ^ = 

Equations (2.31) and (2.32) are solved using Bubnov- Gal erkin 
finite element method (See Appendix A) . 

At the I-th iteration between successive time steps let T^ 
and T^ be the temperature values computed at the interface from 
subdomains l and 2 respectively. Then the guessed gradient 
boundary condition X is updated using Uzawa' s algorithm as 
follows . 






(2.33) 


+ p (T^ - tJ) / L 

As before values of L = 1 and p = 0.2 have been used in this 
study. Figure 2.4 shows a flow chart for implementing Uzawa's 
algorithm. 


2.4.6 Single Cylinder With 9-Subdomains 

The full domain is divided into 9 -subdomains with one of them 
containing the circular cylinder. The central subdomain contains 
the cylinder. All other subdomains are simple rectangles. 

The governing differential equation, initial conditions, 
boundary conditions are given below for the full domain. 

1 2 

Governing differential equation: T^ = ^ V T 

0<x<6, 0<y<6, t>0 (2.34) 

Initial condition: T = 0 at t = 0 

Boundary conditions: x=0,0^y:s6;T = 0 

X - 6, 0 = y = 6; g = 0 
y = 0 and 6, 0ixs6;^=0 

On the cylinder surface, we consider a constant heat flux 


boundary condition 

- 1 , 


dl 

8n 


where, n is a unit outward drawn normal on the cylinder surface. 

The mathematical problem solved in each subdomain is given 
below . 

Subdomain 1 


Governing differential equation: T^ = ^ '7‘‘ 


Initial condition: 
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aT 

an 
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II 
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aT 

an " 

0 

y = 1 / 
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S 

X 

< 

1; 

aT 

an 

x^ 

^2 


(2.35) 


Subdomain 2 




1 2 

Governing differential equation: ^ ^ 

0<x<l, l<y<5, t>0 

Initial condition: T = 0 at t = 0 

Boundairy conditions: x = 0, l£y:£5;T = 0 

X = 1, 1 * y = 5.- 55 = ^3 

tr-i 

y - 1, u - X - 1, an ~ " 2 


aT 


y = 5, 0 = X = 1; 


(2 


Subdomain 3 

Governing differential equation: 
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Pe 


2 

V T 



0<x<l, 5<y 
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6, 
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Initial 

condition: T = 0 

at 
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Boundary 

conditions: x = 0, 

5 

< 

y 

< 

6 ; 

T = 

0 


X = 1, 

5 

< 

y 

< 

6; 

ST 

Sn 



y = 5, 

0 


X 

< 

1; 

ST 

Sn ~ 



y = 6, 

0 

:< 

X 

< 

1; 

Sn 
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Subdomain 4 

Governing differential equation: T 
l<x<5, 5<y<6, 


t 

t > 


Pe 


v^T 


0 


Initial condition: T = 0 at t = 0 

Boundary conditions: x = l, 5:£y^6; 

x=5, 5^y^6; 
y=5, l£xs5; 
y=6, l^x^S; 

Subdomain 5. 

Governing differential equation: = 

5<x<6, 5<y<6, t>0 


ST 

Sn 

ST 

Sn 

ST 

Sn 

£1 

Sn 


= 0 


1 2 
^ V T 
Pe 


(2 


(2 


(2 


.36) 


.37) 


.38) 


.39) 


35 


Initial condition: 


T = 0 at t 


Boundary conditions: x=5, 5^y^6;|^=A 


6, 5 :£ y 


£1 

dn 

dT 


iH = ° 


r- r- ^ dT 

y = 5, 5 ^ X ^ 6; ^ 

dT 
8n 


y=6, 5 ^ ^ 6 ; 


= 0 


Subdomain 6 


1 2 

Governing differential equation: ^ V T 


5<x<6, l<y<5, t>0 

Initial condition: T = 0 at t = 0 


LTJ 

II 

1 

< 

y 


5; 

aT 

an 


X = 6, 

1 

;< 

y 

< 

5; 

aT 

an " 

0 

y = 1 - 

5 

< 

X 

< 

6; 

£1 
an ~ 

0 
H H 

1 

LO 

II 

5 

< 

X 

< 

6; 

£1 

an 

- ^8 


Subdomain 7 


1 2 

Governing differential equation: ^ V T 


5 < X < 6 , 

O 

A 

< 

1 

/ 

t 

> 

0 


Initial condition: 

II 

O 

at 

t 

= 

0 




Boundary conditions : 

X = 5, 

0 

< 

y 


1; 

aT 

an 



X = 6, 

0 


y 

< 

1; 

aT 

an 

0 


o 

II 

5 

;< 

X 

< 

6; 

£r 

an 

0 


y = 1, 

5 


X 

< 

6; 

aT 

an 

^10 


Subdomain ^ 

Governing differential equation: 
1 < X < f 

Initial condition: 


1 2 

T = — V T 
t Pe 


O 

A 

< 

1, 

t 

> 

T = 0 

at 

t = 

0 


X = 1, 

0 

^ y 

:< 

1; 

X = 5, 

0 

^ y 


1; 


dT , I 

5n " ■ 
aT , I 

an " " ^11 


(2.40) 


(2.41) 


(2.42) 
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y = 0, 1:£X:S5; 
y = l, 1 :£X:S 5 ; 


£r 

an 

aT 

an 


12 


Subdomain 9 


. 1 2 

Governing differential equation: T = =— V T 

U i:r0 



1 < X < 5 , 

1 < y < 

5 

f 

t 

> 

0 


Initial 

condition: 

T = 0 at 

t 

= 

0 




Boundary 

conditions : 

X = 1, 1 

< 

y 

< 

5; 

3T 

Sn 

- ^3 



X = 5, 1 

:< 

y 


5; 

3T 

Sn " 

-9^ 



y = 1, 1 

:< 

X 

< 

5; 

ST 

Sn " 

- X^ 
^12 



y = 5, 1 

:< 

X 

< 

5; 

ST 

Sn " 



(2.43) 


On the cylinder surface, 


aT 

an 


= 1 . 


Equations (2.35) to (2.43) are solved using Galerkin finite 
element method. At the I-th iteration between successive time 

steps, let and ’^ 22 ' "^23 '^24' "^34 '^35' ”^45' *^46 

and and T^g, Tgg, Tg^ and Tg^g, and Tg^^, Tg^^ 

and Tg^, Tgg, Tgg, Tgg and temperature values computed at 

the interfaces from subdomains 1, 2, 3, 4, 5, 6, 7, 8 and 9 

respectively. In the notation T^, x represents subdomain number, 
y represents interface number and I represents iteration number. 
The guessed gradient boundary condition X is updated using Uzawa' s 
algorithm as follows 
xi 

J 

x3 

x: 
x; 


I+l 

T 

- 

+ 

P 

(T^ 

^11 

- T^ ) 
■^81^ 

/ L 

I+l 

2 


+ 

P 

(T^ 

^12 

- T^ ) 
22^ 

/ L 

I+l 

3 

= ^3 

+ 

P 

(T^ 

^•^23 

- ) 
93^ 

/ L 

I+l 

4 


+ 

P 

(T^ 

^34 

- T^ ) 
24^ 

/ L 


(2.44) 
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I+l 
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^12 
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P 

(T^ 

^812 
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2.5 VELOCITY FIELD FOR 1 AND 5 CYLINDERS 

In the front problem, the velocity field is prescribed as u = 
1 and V = 0, a constant. For flow past an array of cylinders, the 
flow is assumed to be Darcian, i.e., flow is assumed to be 
governed by Darcy's law. 

u = - ^ VP (2.45) 

The fluid is taken to be incompressible and so V.u = 0. The 
boundary condition u.n = 0 is applied on all impermeable surfaces. 
For flow past a single cylinder these equations can be solved 
using potential theory. This yields, 

u - i V = 1 - ^ , where z=x+iy, i= 1-1 

and u and v are cartesian components of velocity. For flow past 
an array of cylinders the permeability k appearing in Darcy law is 
made a spatial variable with a value of unity in the fluid region 

_7 

and a small value, around 10 over the area occupied by the 
cylinder. The velocity problem is then solved numerically in 
terms of a stream function \ji which is related to velocity as u = 
{rjj , -ip ) . the equation governing ili can be shown to be 

y ^ 

V.i ViJj = 0 (2.46) 

with the approach flow condition ~\p^) = (1, 0) i.e. i{/ = y 

ahead of the cylinder array. 

In the present study calculations have been carried out on a 
square domain in which a single cylinder or array of cylinders is 
present. This is shown in Figure (2.8). With reference to this 


figure the boundary conditions for Equation (2.46) are 

x = 0, }/) = y; y = 0 , ^ = 0 

y = L , Ip = h and X = L, = 0 (2.47) 

Equation (2.46) has been solved by a control volume finite 
difference method with harmonic averaging for k at the interface 
between the cylinder and the porous region. fhe grid size of 101 
X 101 has been found to adequate to obtain grid independent 
velocity profiles. Streamline pattern is shown in Figure (2.9). 

2.6 GRID GENERATION 

Partial differential equation based methods are commonly used 
for structured grid generation in a wide variety of applications. 
These techniques may be termed as elliptic, parabolic or 
hyperbolic, depending upon the characteristics of the grid 
generation equations . In most applications the grid generation 
equation are also transformed on the rectangular computational 
domain and solved along with the governing equations of the 
physical problem. Elliptic grid generation methods offer good 
grid smoothness and spacing control, and have been employed in 
this work. Elements generated by grid generation are 

quadrilaterals. These have been sxibsequently subdivided into 
triangles for finite element calculations. 


u= U 
v= 0 

p= 0 





CHAPTER 3 


TESTING OF COMPUTER PROGRAM 


3.1 NUMERICAL DETAILS 

In the present study, a variety of grids have been used 
depending on whether the problem studied was front movement or 
flow past a cylinder dnd the Peclet number. In the absence of 
cylinders, the grid used is uniform and identical six-noded 
triangular elements are used. Figure 3.1 shows the grid used in 
this present study with 216 six-noded triangular elements for 
problems with Pe 100. The number of nodes used is 481. The 
region size considered here is 18 x 6, dimensionless units. For a 
simulation with one cylinder the grid used is necessarily 
non-uniform (See Figure 3.2). For a region size of 6 x 6, 256 
elements and 576 nodes with 32 nodes on cylinder surface are used 
for Pe i 10 and for Pe 2: 50 region size used is same, 576 elements 
and 1248 nodes with 48 nodes on the cylinder surface are used. 
Owing to the use of isoparametric elements the representation of 
the cylinder is exact. 

The operator splitting technique is explicit in time. The 
maximum permissible time step is found by using the courant number 
restriction, C 1. In the worst case C = 1 where 


C = max 


U At V At 

max max 


(3.1) 


Ax ' Ay ^ 

Here the velocities u, v at all nodes and Ax, Ay for all the 
elements are known before hand. Applying Equation (3.1) to all 
the nodes, the maximum permissible At over the whole mesh can be 
found. Since the definition of C in Equation (3.1) is only 
approximate, a value of 50% of that maximum permissible time step 
has been used for the present calculations . For the front 
movement problem. At = 0.5 for all Peclet numbers. For flow past 
single cylinder. At = 0.095 for Pe ^ 10 and At = 0.05 for Pe > 10 
for the grid shown in Figure 3.2. 
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Figure 3.1 Finite Element Mesh of the Computational Domain 
without Cylinder 
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Figure 3 . 2 Finite Element Mesh of the Cotiputational Domain with 
Single Cylinder 




3.2 FRONT MOVEMENT PROBLEM 

In this section, we consider steady flow between two parallel 
insulated surfaces and a step increase in temperature at the 
inflow plane. The resulting problem is one -dimensional . Operator 
splitting with universal limiters and finite element solutions 
♦over the flow domain are compared with the analytical solution of 
the problem (see Appendix C for the analytical solution) . 

Fronts will develop at high Peclet numbers and will propagate 
in steps of time intervals over the grid spacings. At low Peclet 
numbers the transport is diffusion dominated. Owing to this^ clear 
fronts are not seen in low Peclet number flows . 

The governing differential equation for propagation of a 
thermal front is 

T. + u.VT = ^ V^T (3.2) 

t Pe s 

along with boundary conditions 

0 £ y :£ L 
0 £ y £ L 
0 ^ X :£ L 

where L is non-dimensional length of domain. 

Further, the initial condition is, T = 0 at t = 0 . 

For the analytical solution of Equation (3.2), the outflow 
boundary condition used is, x — > oo, T — > 0. Two time levels t = 5 
and t = 10 are considered for comparing the numerical and 

analytical solutions . 

For Pe = 0.1 and at t = 5 and t = 10 a gradual temperature 
profile moves through the fluid domain (Figure 3.3). The 
analytical solution deviates from operator splitting and finite 
element method after some distance (x = 11) because the analytical 
solution has the different far field boundary condition namely T 
— > 0 at X — > CO. However, the numerical methods compare closely 
with one another. 

At Pe = 1.0 all the 4 methods give identical answers (Figure 
3.4) . At Pe = 10 at both time levels (t = 5 and t = 10) a sharp 
front is clearly seen (Figure 3.5). Here operator splitting and 
analytical methods coincide for all values of x. The conventional 


T = 1, at X 


0 


T = 0, at X = L 
aT 


an 


0 at y = 0, L 















finite element results deviate at positions of sharp gradient 
because of round off errors during inversion of the unsyrametric 
matrix. Hence, the Bubnov-Galerkin method is not directly 
suitabl,e for highly convective flows . 

The same trend as at Pe = 10 is maintained at Pe = 100. The 
performance of operator splitting is very clear compared to that 
at Pe = 10. The fronts at Pe = 100 are shown in Figure (3.6) . 
Operator splitting and operator splitting with limiters yield 
identical answers and the latter has not been explicitly shown. 

3.3 FLOW PAST A SINGLE CYLINDER; COMPARISON WITH BOUNDARY LAYER 

METHOD 

A heated cylinder placed in Darcian flow will develop a 
thermal boundary layer over its surface at high flow rates . 
Boundary layer theory is valid at steady state and for single 
cylinder problems at high Peclet number flows, it can then be used 
to deteirmine heat transfer rates from the cylinder. Bounda 2 ry 
layer theory is not valid during transients, low Peclet numbers 
and in cylinder array problems. Boundary layer solution for flow 
past a cylinder is presented in Appendix D. Using this approach 
we get local Nusselt number over the cylinder surface and the 
average Nusselt number. 

Figure (3.7) shows steady Nusselt number variation on the 
cylinder surface at Pe = 0.1, Pe = 1.0 and Pe = 10. The operator 
splitting and finite element method results are seen to coincide, 
but deviate from boundairy layer theory. At Pe = 10, the boundary 
layer results coincide with those of operator splitting, but 
deviate from the finite element results. In this case, finite 
element method requires a fine mesh for getting -cfr greater 
accuracy. From this we conclude that operator splitting is 
efficient even on coarse mesh at high Peclet numbers. 

Figures (3.8) shows steady Nusselt number variation on the 
cylinder surface at Pe = 50 and Pe = 100. The behaviour of 
operator splitting and operator splitting with universal limiters 
are somewhat different from boundary layer theory. In particular, 
the variation of Nusselt number with angle is non-monotonic when 
numerical methods are used. We can see whether this result is 
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physical or unphysical by using a fine mesh. This is discussed in 
Chapter 5 in the context of domain decomposition. Figures (3.7) 
and (3.8) and Table (3.1) show that operator splitting and 
operator splitting with limiters are generally accurate methods at 
high Peclet number flows, compared to the full finite element 
method (FEM) for a fixed number of elements. 







(too 003 1:20 LOO zsn. 3li4 


Radians 


Figure 3.7 Variation of Steady Nusselt Number 
Surface for Pe = 0.1,1 and 10 


on the Cylinder 
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Table 3.1 Steady Nusselt Number Values for Different Methods 


Method 

Pe = 0.1 

Pe = 1.0 

Pe = 10 

Pe = 100 
19 X 19 

Pe = 100 
23 X 23 

BLM 

0.2645 

0.8368 

2.646 

8.368 

8.368 

FEM 

0.370 

0.922 

3.762 

20.071 

19.844 

OS 

0.370 

0.893 

2.668 

5.382 

5.287 

UL 

0.370 

0.893 

2.668 

5.388 

5.156 



CHAPTER 4 


FLOW PAST 1 AND 5 CYLINDERS; FULL DOMAIN SIMULATION 

In the present chapter, the flow field for a single or an 
array of cylinders is determined using Darcy' s law and is 
discussed in Chapter 2 . There is no separation of flow and there 
is no importance for the hydrodynamic boundary layer. Here we 
concentrate only on the thermal boundary layer and associated heat 
transfer during the transient as well as steady state. 

4.1 HEAT TRANSFER FROM A SINGLE CYLINDER 

For a single cylinder, we compare the results obtained by 
finite element method, operator splitting, operator splitting with 
universal limiters and boundary layer method. For low Peclet 
number the domain size is 6 x 6 units, number of elements is 256 
and number of nodes is 576. For high Peclet numbers the domain 
size is the same the number of elements is 576, the number of 
nodes is 1248 and the mesh is fine. Here we find the mesh size by 
calculating the minimum thermal boundary layer thickness and 
ensure that within this thickness, we have sufficient number of 
elements. Time steps used in present case for Pe ^ 10 is 0.095, 
for Pe i 10 is 0.05 and 0.03. For Pe = 50, and Pe = 100 we used 
two different types of grids (See Figures 4.1 and 4.2) . 

At low Peclet numbers all methods give identical answers at 
steady state except the boundary layer method (Figures 4.3 and 
4.4) . At higher Peclet number flows the classical finite element 
solution deviates from the other methods (see Figures 4.5 and 
4.6) . The boundary layer method is invalid at low Peclet number. 
Its predictions match numerical results at intermediate Peclet 
numbers. At high Peclet numbers the local Nusselt number 
predicted using boundary layer approximation fails to show the dip 
observed with numerical schemes. This issue is discussed again, 
later. The average Nusselt number at small time is large and 
decreases rapidly to reach steady state. For a constant heat flux 
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Figure 4 . 
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Finite Element Mesh with Single Cylinder (type 1) 
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2 Finite Element Mesh with Single Cylinder (type 2) 
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Figure 4.3 Variation of Steady Nusselt Number on the Cylinder 
Surface for Pe = 0.1, 1 and 10 
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Figure 4 . 6 Variation of Average Nusselt Number with time for 
Pe = 50 and 100 





specified over the cylinder surface, = 1/T, where T is the 
local wall temperature . Initially the average temperature of the 
cylinder is small and increases gradually with time. Time shown 
in plots for average Nusselt number is defined as, t = at/L^ and 
is the non-dimensionalized time scale. 

At low Peclet numbers the steady average Nusselt number is 
low compared to its value at high Peclet number. This is because 
at high Peclet number convection dominates transport . A greater 
amount of heat is lost from the cylinder and the temperature of 
cylinder decreases. At low Peclet number (Pe ^ 10) the finite 
element method, operator splitting and operator splitting with 
universal limiters coincide (Figure 4.4). At high Peclet numbers 
(10, 50, 100), the finite element method deviates from operator 
splitting and operator splitting with universal limiters (See 
Figure 4.6). For completeness, the no flow problem, Pe = 0.0 is 
also presented here (Figure 4.7) . There is no steady state in the 
absence of flow and Nusselt number goes to zero at long time. 


4.2 HEAT TRANSFER FROM FIVE CYLINDERS 

In case of 5 cylinders, a boundary layer solution is not 
possible to account for the inter- cylinder influence. Instead, we 
have compared the three methods namely finite element method, 
operator splitting and operator splitting with universal limiters. 
The domain size is 11 x 11 units. The number of elements is 808 
and the number of nodes is 1780. The five cylinders are located 
at (5.5, 5.5), (2.5, 2.5), (2. 5, 8. 5), (8. 5, 2. 5), (8. 5, 8. 5), as 
seen in Figure 4.8. The maximum allowable time step is calculated 
using the Courant number condition. The time step is calculated 
as At = 0.09 for all the three Peclet numbers namely, ,0.1, 1 and 
10 . 

The thrust of the calculation is to determine the effect of 
the 4 surrounding cylinders on the central cylinder. At Peclet 
numbers of 0.1 and l the steady Nusselt number value is less than 
that a single cylinder but at Pe = 10 the difference is small (See 
Figure 4.9) . The reason for this behaviour is the following at 
low Peclet number, diffusion is main mechanism of transport of 
energy. All 5-cylinders have a constant heat flux boundary 
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condition. Due to the other cylinders located around the central 
cylinder, the temperature gradient decreases and the heat transfer 
from central cylinder reduces. This in turn affects the steady 
Nusselt number of this arrangement. 

From Figure (4.10), variation of steady Nusselt number with 
angle is comparable upto Pe = 10 with all three methods. 
Simulation with higher Pe was not possible owing to memory 
limitation on the HP system. 




Figure 4.8 Finite Element Mesh of the Computational Domain with 
Five Cylinders 
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igure 4.9 Variation of Average Nusselt Number with time for 
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Figure 4.10 Variation of Steady Nusselt Nuniber 
Surface for Pe =0.1, 1 and 10 for 
Five Cylinders 


on the Cylinder 
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CHAPTER 5 


RESULTS OBTAINED USING DOMAIN DECOMPOSITION 


5.1 MOVEMENT OF FRONT PROBLEM 

In the present study, the domain divided into 2, 4 and 9 

subdomains lias been considered. The formulation of domain 
decomposition is discussed in Chapter 2. The full domain has the 
size of 18 units- length, 6 units width. The time step used is 0.5 
for all Peclet numbers. With 2-subdomains we treat only one 

interface and so it is faster than the computations for 4 and 9 
subdomains . In this present case results have been compared with 
the analytical solution. Agreement between analytical and domain 
decomposition is excellent (See Figures 5.1, 5. 2, 5.3, 5.4, 5.5 

and 5.6) . The efficiency of domain decomposition calculation also 
depends upon the type of division. In case of 2 - subdomains , the 
partition is meaningful since the transport of energy from left to 
right. In case of 4 and 9, subdomains the domain is divided 

horizontally as well without a meaningful physical mechanism being 
present in this direction. Convergence rates are fast only in case 
of physically meaningful division of the domain. These are clear 
from Table 5.1, where CPU time is given. The CPU times for domain 
decomposition is more compared to full domain calculations. If 
these calculations were to be implemented on a parallel computer, 
the CPU time is expected to fall below that for a full domain. 

5.2 FLOW PAST A SINGLE CYLINDER, EFFECT OF OUTFLOW PLANE 

5T 

The outflow boundary condition used in the present work is ^ 
= 0, where n is the outward drawn normal at that plane. This 
condition means that, temperature from that plane onwards is 
constant along the streamline direction. This boundary condition 
is valid if the outflow plane is faraway from the cylinder 
surface. Solving the governing equation over large domains on 
0 xisting computers is a problem. Domain decomposition is a 
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variation of Temperature along Flow Directxon wrth TVo 
subdomains for Pe = 0.1 an 


Figure 5 . 1 





































Figxire 5 . 5 Variation of Temperature along Flow Direction with. 
Nine Subdomains for Pe = 0.1 and 1 











convenient approach available at this point to circumvent this 
difficulty. 

Originally a domain of size 6x6 was used for a cylinder 
diameter of 2 -units. Compared to this diameter a length of 6 
units is not large enough. This results in spurious reflection of 
errors from the outflow plane. At low Peclet numbers the 

transport of energy is due to diffusion and there is no effect of 

the outflow plane. In case of high Peclet number flows the effect 
of outflow plane is greater, because transport of energy is 
predominantly due to convection. 

From Table 5.2 the effect of the location of the outflow 
plane is clear enough. 

For the two domain problem, the domain size is 12 x 6. For 
low Peclet numbers, 544 elements and 1276 nodes and for high 
Peclet number 1224 elements and 2580 nodes have been used. Figure 
(5.7) shows Nusselt number as a function of angle e for a domain 
of size 6x6 and another for a domain of size 12 x 6 units with 

2 -subdomains each having 6x6 units of region size. For Pe < 10 

both curves coincide showing that the outflow location is not 
critical. For Pe > 10, the curves deviate from each other but the 
overall variation is the same. From this figure, the outflow 
effect is clear, because there is a change in steady Nusselt 
numbers over the cylinder surface. Figure 5.8 shows the 
variations of Nusselt number with time. The isotherms patterns 
are shown in Figure 5.9. From this figure, it can be observed 
that at low Peclet numbers the outflow condition is valid and at 
high Peclet numbers the outflow condition is not exactly 
satisfied. 

5.3 9 - SUBDOMAIN SIMULATION FOR FLOW PAST A CYLINDER 

In this section, we present results for a domain divided into 
9 sub -domains with one cylinder at the centre of the domain (See 
Figure 5.10) . The formulation of domain decomposition is discussed 
in Chapter 2 . The full domain has the size of 6 x 6 units . The 
number of elements is 256 and number of nodes is 576 for Pe = 0.1, 

1 and 10. The time step used is O.l. Here we compare the 
variation of steady Nusselt number with angle 6 with that of the 
full domain. At low Peclet numbers the variation of steady 
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Table 5.1 CPU Times in Minutes for Interface Error = l.Oe-04 : 
Without Cylinders 


Peclet 

Numbers 

Full domain 

2 -Subdomain 

4 - Subdomain 

9 - Subdomain 

0.1 

1.828 

1.18 

1.26 

1.03 

1.0 

1.74 

1.85 

3.20 

2.98 

10 

2.35 

3.33 

9.66 

9.00 

100 

2.28 

10.91 

42.21 

43.66 

Table 5 . 2 

CPU Times in Minutes for 
Without Cylinders 

Interface Error = l.Oe-06 

Peclet 

Numbers 

Full domain 

2 -Subdomain 

4 -Subdomain 

9-STibdomain 


0.1 

1.828 

1.51 

4.25 

- 

1.0 

1.74 

3.85 

10.85 

8.81 

10 

2.35 

10.01 

26.28 

26.7 

100 

2.28 

32.70 

113.91 

115.9 


Table 5.3 Steady Nusselt Number Values for Different Finite 
Element Meshes 


Method 

Pe = 0.1 

Pe = 1.0 

Pe = 10 

Pe = 100 

6x6 

0.370 

0.893 

2.668 

5.382 

12 X 12 

0.368 

0.896 

2.662 

5.544 
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Nussslt nuttibsirs with sng'l© Q (S©© Fi^uirs 5.11) is ici©nticsl to 
that of full domain. Th© variation of averag© Nuss©lt number with 
time level is compareti in Figure 5.12. Th© comparison between 
full domain and 9 subdomain calculations is good. This is a good 
check on the domain decomposition. Hence it can be implemented for 
5-cylinders, for large domains and for high Peclet number flows. 

We present an explanation for drop in Nusselt Number around e 
= 1^0° for Peclet numbers greater than 50. The local Nusselt 
number on the surface of the cylinder is subjected to the 
conflicting influences of increase in boundary layer thickness in 
the flow direction and increase in velocity between e = |eo° to 
The increase in the boundary layer thickness reduces the 
local Nusselt number in contrast to this, the increase in local 
velocity increases the local Nusselt number. The breakthrough 
between these two factors will occur at an angle that depends on 
the Peclet number. For Pe s 10 the effect of boundary layer 

growth is dominant. For Pe ^ 50 both effects are simultaneously 
present under non-monotonic variation in the local Nusselt number 
is seen. On the rear of the cylinder {0°s |e| :£ 90°) . The flow 
is decelerated and Nusselt number drops at fairly rapid rate. 










t 


Figure 5 . 12 Variation of Average Nusselt Nuitiber with time for 
Pe = 0.1, 1 and 10 : Nine Subdomains 







r iguic o.A *V3.i i3.tion of Stcs-ciy Nussclt Nunibor on the Cylinder Surface 
for Pe = 100 


Figure 5.1) Variation of Average Nusselt Number with 


time for Pe = 100 





CHAPTER 6 


COMMENTS ON PARALLELIZATION 

To take advantage of the inherent flexibility of the finite 
element method in solving for flows within complex geometry, it is 
necessary to produce an efficient implementation of the method. 
Functional splitting and domain decomposition are two methods 
which will make finite element method very efficient. Both 
functional splitting and domain decomposition are useful for 
parallelization of any problem. 

Two types of parallel computers, those are based on 
architecture are widely used. Those are single instruction, 
multiple-data stream (SIMD) and multiple instruction, 
multiple-data stream (MIMD) . The single instruction multiple-data 
stream computers are extensively used in scientific applications. 
Developing a computer code for single-instruction, multiple-data 
stream computer is difficult and in some cases it is not possible, 
for example when subdomains are not equal in size or shape. 
Another type of parallel machine is a combination of different 
workstations through some software linkage, for example Parallel 
Virtual Machines. 

The important factors to be considered in parallelization 

are : 

1 . Load balancing 

I 

2 . Communication efficiency 

3 , Interface treatment 

Equal load on all processors will give good performance on a 
parallel computer. This is a difficult task and maintaining this 
is cumbersome in some cases because of unequal domains . 
Communication efficiency depends on the number of interface nodes. 
Hence, interface treatment which depends upon the algorithm used 
for interface treatment is crucial for improving parallelization 

efficiency. 
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CHAPTER 7 


CONCLUSIONS AND SCOPE FOR FUTURE WORK 

7.1 CONCLUSIONS 

Ihe performance of operator splitting algorithm was tested in 
simple as well as complex geometries against the Bubnov-Galerkin 
method. Results were compared against the analytical solution of 
problems which were available in some cases . It was found that 
for low Peclet number flows, the performance of operator- splitting 
was identical to performance of other methods. For high Peclet 
numbf'r caso:^ only the results by operator splitting came close to 
the analytical solution. Universal limiters were implemented for 
removing interpolation oscillations arising in high Peclet number '■ 
flows. Operator splitting with universal limiters is 
relatively free from the upwinding errors and successfully 
simulates flows at high Peclet numbers as well as in complex 
geometries . 

The performance of domain decomposition was tested in simple 
as well as complex geometries against the full domain problem. 
Results were compared with analytical solution which were 
available. The performance of domain decomposition depends mostly 
on the orientation of the interfaces, type of physical problem and 
interface treatment. This is clear from 2-domain, 4-domain, 

9 -domain problems with no cylinder and a single cylinder. The 
speed of domain decomposition depends on interface treatment and 
convergence criterion. For a convergence less than 10 , there is 

no Lmprovoment in the results and only the number of iterations 
increases. With the help of domain decomposition we have tested 
the effect of outflow boundary on heat transfer from a cylinder 
placed in flow. 

7.2 SCOPE OF FUTURE WORK 

The following issues must be addressed in the future. 

1. The performance of operator splitting and domain 
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/ 27 


decomposition for Navier-Stokes equations must be studied. 

2. Uzawa's algorithm should be compared with conjugate gradient 
method for interface treatment to check for possible 
improvement . 

3 . The computer code developed in the present study must be 
implemented on a parallel machine to study improvements in 
CPU time. 

4. Sensitivity of the solution to various meshes must be 
examined . 

Simulations at very high Peclet numbers must be carried out. 


5 . 
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APPENDIX A 


FINITE ELEMENT METHOD 

The finite element method is a powerful computational 
technique for the solution of differential and integral equations 
that arise in various fields of engineering and applied science. 
The finite element method is particularly advantageous for complex 
geometries and gradient boundary conditions especially when 
compared to the finite difference method. The major steps in 
finite element method applied to a typical problem are as follows: 

1. Discretization of the domain into a set of finite elements 
(mesh generation) . 

2 . Weighted - integral or weak formulation of the differential 
equation to be analysed. 

3 . Development of the finite element model of the problem using 
its weighted integral or weak form. 

4. Assembly of finite elements to obtain the global system of 
algebraic equations . 

5 . Imposition of boundary conditions . 

6. Simultaneous solution of the algebraic equations. 

7. Post-computation of solution and quantities of interest. 

Steps 

1 . Mesh Generation 

In the present work, six-noded isoparametric triangular 
elements are used for accommodating the cylinder configuration 
correctly. Elliptic grid generation procedure is used for 
smoothening grids . 

2 . Weighted Integral or Weak Formulation 

The governing differential equation is 

T^ + u.VT = V^T (A.l) 
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where is time derivative 

„ . ^ a t a 

V IS 1 T— + 1 ^ 

au ay 

2 2 

„2 . a^ a^ 

V IS — ^ + =r 


ax ay^ 

The symbols T, u are temperature and velocity respectively. 

To develop the weak form, we multiply Equation (A.l) by weight 
function w and integrate the equation over the element doraain^^ /le 
Hence, 


S 

ne 


[T^ + U.\7T 


V^T] W dX = 0 


(A. 2) 


where dX = fdx dy ■ 

After simplification of Equation (A. 2) we get 

S [T.W + (U.VT) W - (VW.yT)] dX 
ne ^ 

-fe fl^Wds = 0 

Pe an 


(A. 3) 

0 at t = 0 


Equation (A. 3) subjected to initial condition T 

9T 

and boundary conditions are ^ + hT = q on Fe where £2e is the 

element domain and Fe is the boundary of that domain. 

4 . The dependent variable T is approximated over a typical 

element n® by the expression 

n 

T (x,y) « T'" (x,y) = X T^'^'Ct) N^‘^(x,y) (A. 4) 


T® (x,y) = X T.®{t) N.®(x,y) 
j=l ^ ^ 


where j = 1, 


n, n is number of nodes per element. Here N. 


is shape function and is function of x and y only, and Tj is 
function of time. T^® denotes the values of the function T®(t,x,y) 
at selected number of points called nodes . 

After imposing boundary conditions in Equation (A. 3), we get 

1 


S [T.W + (U.VT) W - ^ (VW.VT)] dX 
Qe ^ 


Pe 


f h T W d s 
Fe 


1 

Pe 


f q W d s = 0 
Fe 


(A. 5) 


Substituting Equation (A. 4) in (A. 5) we get, 

n n 

S [T^ W + (U.V X Tj® Nj®) W - ^ (VW.V X 


Qe 


j=l 

n 


j=l 


T.® N.®)] dX 
3 3 


+ ^ f h X T.® N.® Wds - f q W ds = 0 
^® Fe j=l ^ ^ 


Fe 


(A. 6) 





This equation must hold for any weight function W. in 
Bubnov-Galerkin method weight functions are equal to shape 
functions i.e. W = . Since we need n independent algebraic 

equations to solve for the n unknowns, T. T We 

-L /L jQ. 

choose n independent functions for W namely W = For 

each W, we obtain an algebraic equation among We 

label the algebraic equation resulting from substitution of 
for W into Equation (A. 6) as the first algebraic equation. The 
i-th algebraic equation is obtained by substituting W = into 

Equation (A.6) . 


Equation (A.6) becomes, 

Z (M?. T® + K®. T®) = Q® 

13 J 1: 1 1 

where the coefficients M? . , K?. and Q® are 

13 ^3 ^ 

Mf.=f N® N® dX 
, ^ fie ^ ^ 


(A. 7) 


K?. = S [(U.VN.®) N.® - I- (VN.®. VN.®)] dX 

-LJ fig J i ir-e X J 


+ f h N. ® N ® ds 

^® re ^ J 


Q.® _ 1 

1 Pe 


f q N. ds 

re 


In matrix notation. Equation (A. 7) takes the form 
[M®] [T®] + [K®] [T®] = [Q®] 


(A. 8) 


0 0 
where the matrix [M ] is called as mass matrix and [K ] is called 

as stiffness matrix or conductivity matrix. 

Implicit# finite difference scheme has been used in the 
present work ' for temporal approximation. After assembly of all 
these element matrices, we get the global matrix. This can be 
inverted by using direct methods or indirect methods. A direct 
inverter using Gaussian elimination with full pivoting has been 
used in the present work. 

The integrals appearing in the definitions of the matrix 
entries are computed using Gaussian-quadrature. We have used 7 
Gauss points in six noded isoparametric elements. These are, 


X : 0, 1, 0, 0, 0.5, 0.5, 0.3333 
3 



y^: 0, 0, 1, 0.5, 0, 0.5, 0.3333 

The associated weights used in present work are, 

1 1 1 11 1 
40 ' 40 ' 40 ' 15 ' 15 ' 15 ' 4.444 

The quadratic shape functions defined over a triangle are listed 

below. 


Nl 

II 

H 

H 

- 1) 

^2 

= (2L2 

- 1) 

^3 

= L3 (2L3 

- 1) 

^4 

= 4 L 1 L 3 


^5 

= 4L1L2 


^6 

= ^^2^ 



where , 


= 1 - e - 11 

s = « 

Lj = 7) 

where t) are natural coordinates on a unit triangle. 



APPENDIX B 


METHOD OF CHARACTERISTICS FOR ADVECTION STEP IN 
OPERATOR-SPLITTING ALGORITHM 


The advection step in the OS algorithm generates is a 
hyperbolic equation. Consider the advection step of Equation 
( 2 . 1 ) , 


+ u. VT = 0 

In full form the above equation is 


(B.l) 


+ u + V = 0 

The above equation is simplified using streamline 
coordinates. We choose a curvilinear coordinate system (tp, -<P) 
locally at each node such that V is aligned with the net 

composite velocity at that node, i.e. ip is the streamline passing 
through the node [See Figure B.l] . Functions ip and <p are 
orthogonal to each other. Then, it can be shown that 


where U 


= I 


and 


u * V Ty = U 

T - 21 
rp dip 

Equation (B.l) now reduces to. 


u 


+ V 


Tt + U T^ = 0 , 

The solution of this equation is 


(B.2) 


T ( Ut- ip ) = constant (B.3) 

Since the grid surrounding each node is small, we assume the 
streamline ip to be locally straight over an element. 

The solution of the corrector step requires that the 
temperature at 'Q' at n-th time step be transferred to the point 
'P' at the (n+l)-th time step (see Figure B.l) i.e., 

Tp"^^ ( U(t + At) - ip) =Tq { TJt - { ip - tip)} (B.4) 
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where = U.At 

If the advection step is solved after the diffusion step, 
refers to the solution obtained at the end of calculation of 
diffusion step. Here for finding the temperature value at 'Q' we 
use second order Lagrangian interpolation. The local master 
element coordinates ri are found from global x and y coordinates 
by using isoparametric relationship between v and x, y. These 
isoparametric relations are used for mapping the original domain 
in the physical plane to a unit triangle. These relationships are 


n 



X = Z X. N. 
j=l 3 D 

(?, T?) 

(B.5) 

n 



y = Z y . N . 
j=l 

V) 

(B.6) 


where n is number of nodes per element, are shape functions and 
are functions of ^ and V- Xj and y^ are coordinates of nodal 
points. Equations (B.5) and (B.6) are solved using Newton-Raphson 
method (See Appendix E) to determine ^ and v as solutions of 
Equations (B.5) and (B.6). 



Figure B.l Streamline Passing Through the Node. 



APPENDIX C 


ANALYTICAL SOLUTION 


ANALYTICAL SOLUTION FOR EQUATION (2.1) 

The advection-dif fusion equation in one dimensional form with 
unit velocity is given by 


T.. +T = ^ T 
t X Pe XX 


= 0 , 


(C.l) 

T = 0 and the boundary 


subjected to initial condition: t 
conditions : 

X = 0, T = 1 
X — > 00, T = 0 

This represents a physical problem in which a hot fluid with 
unit velocity is introduced suddenly (at t = 0) at the inlet of a 
flow domain. The fluid in the physical domain (0 < x < oo) is 
initially at a constant cold temperature. If the Peclet number is 
large, a thermal front whose thickness is small will move through 
the flow region with unit velocity. As Pe — > 0, the front will be 
diffuse and energy transport is due to the molecular action of 
thermal diffusivity. 

The closed-form solution of the above equation is, 

2 


T(x,t) = 




00 

s 


e^® [ I (x^Pe/16 t?^)] d7? 


I 


4t 

Pe 


where t? = 4t/Pe . 

The above integral is well behaved and is evaluated using 
Simpson's rule for numerical integration using 301 points. 
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APPENDIX D 


BOUNDARY LAYER SOLUTION 


steady flow past a single cylinder at high Peclet numbers (Pe 
» 1) will lead to the formation of boundary layer on its surface. 
This is shown in Figure (D. 1) . All changes in temperature are 
confined to the region 1 < r < 1+5 and is zero outside this 
region. "^e boundary layer thickness shows on the cylinder 
surface and is a function of ' 6' . If the boundary layer is thin 
i.e. 5/R « 1, where 5 is the theimal boundairy layer thickness, the 
governing equation for T can be simplified. If 5/R « l the 

derivatives parallel to the cylinder wall are small in comparison 

2 2 

B 3 B^ B^ 

to those normal to it i.e. oh ^ ^ . The reduced 

de^ Sr^ 

governing equations are. 


(uT)q = pg ( (D.l) 

where r and 6 are polar coordinates measured from the centre of 
the cylinder and u is now the tangential conponent of velocity. 
From potential theory the tangential conponent of velocity on the 
cylinder surface is u = - 2 sin B for flow past a single cylinder. 

Equation (D.l) is solved by first assuming a profile for T in 
terms of r with 5 as a parameter. Boundary layer thickness 5 is 
then determined from an integrated form of Equation (D.l) . 

The form for T assumed in the present work is 


T = 


( 1 - ( 


r-1 

5 


)) 


It satisfies the boundary conditions r = 1, -T^ = l and r = 


1+5, T = 0 automatically. 
5/2. “^he 


The value of T at cylinder surface is 
'5' is obtained by integrating 


equation governing 
1+5 

Equation (D.l) as X dr. After integrating we get the following 

1 

first order ordinary differential equation. 


de 

where <f> 


+ (p cot 0 = 

.2 


■3 cosec 6 
Pe 


(D.2) 
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Equation (D.2) is solved by a fourth order Runge-Kutta 
scheme. Equation (D.2) exhibits singularity at © = 0 and 5 and T 
are unbounded at this location. This point is excluded while 
determining the average Nusselt number of the cylinder by 
Simpson's rule. 



Figure D.l Thermal Boundary Layer on the Surface of the Cylinder. 



APPENDIX E 


NEWTON-RAPHSON SCHEME FOR DETERMINING (?, -n) FROM (x,y) 


WITHIN A TRIANGULAR ELEMENT 


The advection step requires transfer of temperature data from 
an interior point Q of one element to one of its nodes P (Appendix 
B) . Figure (E.l) shows the isoparametric 9-noded Lagrangian 
master element. For finding (^, v) of point Q, the pair of 
simultaneous equations, 


9 

^ ;] =1 

must be solved. 


N. = f^ (e, TJ) = 0 


N. = g 
J ^ 


Q 


(?, T?) = 0 


(E.l) 


(E.2) 


According to Newton-Raphson scheme, we have to 
compute the corrections and to an initial guess using the 
following expressions. 


A?^ = 


At}*^ = 


/«Q Q Q «Q\ 

5 W 


det 


(g^ f^ 




det 


where det = ( g^ f^- f^ g^ ) 


.Q _ 3g 


,Q 




and 


fQ ^ af^ 


ac 


Here we know the x, 


Q' ^Q' 


X . and 


We assume (Cf "n) 
Q 


coordinates at point Q and calculate A^ and Arj^. If the 
calculated values of A^^ and Atj^ are small, then the guessed (^, 
T)) coordinates are correct, otherwise we modify C V using 

following expressions. 


Q 

new 

Q 

new 


‘old 


+ A?^ 


Q A Q 

”013 


93 



This procedure is repeated till convergence is reached. 

The shape functions used for locating pointy Q are associated 
with a 9-noded Lagrangian element and are given as follows. 



= 

1 

- 

1) 

(T? - 1) 

? 7J 

^2 
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) (t? - 1) 

7? 

^3 
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(t? + 1) 

e 7? 

^4 
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1) 
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^5 
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^7 
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Figure E.l 9 Noded Isoparametric Lagrangian Master Element. 




